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GLOBAL ARTIFICIAL BOUNDARY CONDITIONS FOR COMPUTATION OF 
EXTERNAL FLOW PROBLEMS WITH PROPULSIVE JETS* 


SEMYON TSYNKOVt, SAUL ABARBANEL* , JAN NORDSTROM^, VIKTOR RYABEN’KII^ AND VEER VATSAll 

Abstract. We propose new global artificial boundary conditions (ABC’s) for computation of flows 
with propulsive jets. The algorithm is based on application of the difference potentials method (DPM). 
Previously, similar boundary conditions have been implemented for calculation of external compressible 
viscous flows around finite bodies. The proposed modification substantially extends the applicability range 
of the DPM-based algorithm. In the paper, we present the general formulation of the problem, describe our 
numerical methodology, and discuss the corresponding computational results. The particular configuration 
that we analyze is a slender three-dimensional body with boat-tail geometry and supersonic jet exhaust in 
a subsonic external flow under zero angle of attack. Similarly to the results obtained earlier for the flows 
around airfoils and wings, current results for the jet flow case corroborate the superiority of the DPM-based 
ABC’s over standard local methodologies from the standpoints of accuracy, overall numerical performance, 
and robustness. 

Key words, external flow problems, jet exhaust, artificial boundary conditions, difference partials 
method 

Subject classification. Applied and Numerical Mathematics 

1. Introduction. Many typical problems in aerodynamics including those that present immediate 
practical interest, e.g., flows around aircraft, are formulated on infinite domains. It is, however, obvious, 
that any discretization used for solving such problems on the computer must be finite. Therefore, any 
numerical solution methodology for these problems has to be supplemented (or, rather, preceded) by a 
special technique that helps create such finite discretizations. 

A widely used approach to this problem is based on truncating the original flow domain prior to the 
actual discretization and numerical solution. Subsequently, one can construct a finite discretization on the 
new bounded computational domain using one of the standard techniques: finite differences, finite elements, 
or other. However, both the continuous problem on the truncated domain and its discrete counterpart 
will be subdefinite unless supplemented by the appropriate closing procedure at the external computational 
boundary. This is done by using artificial boundary conditions (ABC’s); the word “artificial” emphasizing 
here that these boundary conditions are necessitated by numerics and do not come from the original physical 
formulation. 

The ideal or, in other words, exact, ABC’s are obviously those that would drive the error associated 
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with domain truncation to zero. However, numerically efficient procedures of this kind cannot be attained 
routinely except in model (mostly one- dimensional) problems and therefore, for typical applications one uses 
primarily different approximate rather than exact methodologies. 

The nature of the difficulties associated with constructing he exact ABC's is that in most cases such 
boundary conditions appear nonlocal (in space and also in time for unsteady problems). Although the cor- 
responding computational algorithms are robust and highly accurate, they can be cumbersome and typically 
apply only to rather simple geometries. On the other hand, bhe alternative local approaches that yield 
inexpensive and geometrically universal numerical procedures may often lack accuracy in computations, 
which, in turn, necessitates choosing excessively large computational domains. Basically, higher accuracy 
due to boundary conditions implies that more of the nonlocal nature of exact ABC's has to be taken into 
consideration. As a consequence, to avoid extra complexity due to the nonlocality of boundary conditions, 
most of the modern production algorithms in CFD still employ local ABC’s that are based on simplified 
flow models. The possibility to use local ABC’s comes, as mentioned, at the expense of running the cases 
on large domains. 

Generally, it has been demonstrated theoretically and computationally in both CFD and other areas of 
scientific computing that the treatment of ABC’s may have a profound impact on the overall performance 
of numerical algorithms and interpretation of the results. The literature on various ABC’s techniques is 
extensive, a detailed review can be found in work by Givoli[l, 2], as well as in a more recent paper by 
Tsynkov[3], 

The construction of ABC’s based on the difference potentials method (DPM) by Ryaben’kii [4, 5, 6], was 
an attempt to combine in one technique the advantages relevant to both local and global methodologies, see 
Refs. [7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17]. These boundary conditions employ finite- difference counterparts 
to Calderon’s pseudodifferential boundary projection operators and generalized potentials that have been first 
proposed in work by Calderon [18] and then also studied by Seeley [19]. The DPM-based ABC’s have been 
successfully implemented along with NAS A- developed multigrid Navier-Stokes solvers for the calculation of 
two-dimensional (solver FLOMG by Swanson and Turkel [20, 21, 22]) and three-dimensional (solver TLNS3D by 
Vatsa, et al. [23, 24]) compressible viscous flows around airfoils (NACA0012, RAE2822) and wings (ONERA 
M6). 

In many numerical tests the DPM-based boundary conditions have consistently outperformed the stan- 
dard local methods from the standpoints accuracy, multigrid convergence rate, and overall robustness (they 
allow for a substantial reduction of the domain size while preserving the accuracy and may also speed up the 
convergence of multigrid iterations by up to a factor of three, i.e., they would require only about one third 
of the original number of multigrid cycles for reducing the initial residual by a prescribed factor). Note, the 
standard local boundary conditions for external flows that are r< ferred to above are typically based on one- 
dimensional characteristics analysis for the front or inflow part of the artificial boundary and specification 
of the free-stream pressure and extrapolation of all other quant ities on the rear or downstream portion of 
the outer boundary; this treatment may or may not be supplemc nted by the point- vortex correction [25] for 
the two-dimensional case; an example of geometry in three dimsnsions is shown on Figure 2.1 in the next 
section. 

All the problems analyzed previously in the DPM framework (see the aforementioned references) can 
actually be characterized as “pure” external flows. In this paper , we for the first time incorporate a new and 
essentially different physical demerit into the formulation of the problem; namely , we will consider external 
flows around the configurations with jet exhaust. The problems of this kind have never been studied by means 
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of the DPM before and including this new flow phenomena into the range of admissible formulations for the 
DPM- based methodology substantially enlarges the scope of its capabilities. Moreover, as different flows with 
jets are frequently encountered in aerospace applications, the possibility of computing external aerodynamics 
more efficiently with jet exhaust phenomena taken into account is important for both configuration analysis 
and design. 

The material in the paper is prepared as follows. In the next section we outline the basic DPM-based 
procedure as developed for pure external flows; in the section that follows we describe the changes that are 
necessary for incorporating the jet exhaust flows; then, wc present the numerical results and conclusions. 

2. DPM-based ABC’s: Basic Algorithm. In this section, we essentially reproduce the correspond- 
ing derivation from Ref. [15]. The paper by Tsynkov [16] contains a substantially more detailed account on 
how to set the three-dimensional DPM-based ABC’s. 




y 
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Fig. 2.1. Schematic geometric setup for the ONERA M6 wing; wing on the left is enlarged. 

We consider a steady-state flow of a viscous, compressible, perfect gas past a finite three-dimensional 
configuration. The flow is uniform and subsonic at infinity; it is also symmetric with respect to the Cartesian 
plane 2 = 0. The hydrodynamics equations are discretized and integrated on a grid generated around the 
immersed body(ies). The grid actually defines a bounded computational domain; the ABC’s that would 
close the truncated problem should be set at the external coordinate surface of the grid. Let us denote 
this surface F; for a one-block curvilinear C-O type boundary- fitted grid around the ONERA M6 wing the 
schematic geometric setup is shown in Figure 2.1. 
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The outermost coordinate surface of the grid is designated Ti (see Figure 2.1); it represents the ghost 
nodes (or ghost cells for the finite- volume formulation). Clearl}, when the stencil of the scheme used inside 
the computational domain is applied to any node from T, it generally requires some ghost cell data. Unless 
the required data are provided, the finite- difference system solved inside the computational domain appears 
subdefinite, i.e., it has less equations than unknowns. Therefore, in practical framework the closure of the 
discretized truncated problem means specification of the solution values at the ghost cells. This will be 
done by means of the DPM-based ABC’s so that the boundary data used for the closure admit an exterior 
complement that solves the problem outside the computational domain. As soon as the data in the ghost 
cells have been obtained as functions of the data in the interior ce lls (Fi as a function of F), the corresponding 
relations can be incorporated into the actual solver used inside the computational domain. If, for example, 
this is an iterative solver (very often the case), then one has tc update the ghost cells at each iteration to 
advance to the next “time” step. 

To construct the boundary conditions, we first assume that the flow perturbations against constant free- 
stream background are small in the far field and consider the linearized problem outside the computational 
domain (i.e., outside T). It is important to emphasize that the possibility of far-field linearization (i.e., the 
possibility to retain only the first-order terms with respect to perturbations in the governing equations) 
requires special justification, in particular, when analyzing the transonic flows. We do not present the 
corresponding argument here; a simple asymptotic analysis in the framework of the full potential model that 
justifies the far-field linearization in three dimensions can be found in our previous work [15, 16]. Of course, 
even if we know that the far field is linear, we still cannot say a priori whether or not the linearization 
outside T is possible for a particular configuration of the domains. Clearly, for a very large computational 
domain one can linearize the flow outside T, and as we approach the source of perturbations (the immersed 
configuration), the validity of linearization is verified a posteriori (see, e.g., Refs. [8, 9, 12, 15, 16]). 

We will be considering the entire problem in the framework of the thin-layer equations (as opposed to 
the full Navier-Stokes equations). This simplified flow model still retains all the essential properties pertinent 
to the class of problems that we are studying and at the same time it is less expensive numerically. In the 
far field, we write down the linearized dimensionless thin- layer equations as follows 


(2.1a) 


dp du dv dw 

dx ^ dx dy dz 

du dp 1 / d 2 u d 2 u\ 

dx dx Re \ch/ 2 dz 2 ) 

dv dp 1 /4 d 2 v d 2 v 1 d 2 w \ 

dx dy Re \ 3 dy 2 dz 2 3 dydz ) 

dw dp 1 / 4 d 2 w d 2 w 1 d 2 v \ 

dx dz Re \ 3 dz 2 ^ dy 2 3 dydz ) 

dp 1 dp k (d 2 p <9 2 p\ 1 (d 2 p 3 2 p\l 

dx Mq dx Re Pr \<9y 2 + dz 2 ) f \9t/ 2 + dz 2 ) \ 


0 

0 

0 

0 

0 


where p, u, u, w, and p are the perturbations of density, Cartesian velocity components, and pressure, 
respectively, Re is the Reynolds number, Pr is the Prandtl number, and k is the ratio of specific heats. The 
full dimensionless quantities at infinity are: p 0 = 1, u 0 = 1, v 0 = 0, w 0 = 0, p 0 = System (2.1a) is 

supplemented by the homogeneous boundary condition at infinity: 
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(2.1b) 


u = (p, u, v , w , p) — * (0, 0, 0, 0, 0) 
as (x 2 + t/ 2 + z 2 ) — > oo, 

which corresponds to the free stream limit of the solution. As has been mentioned, the DPM-based ABC’s 
will close the discretized truncated system by providing the missing external boundary data. These data will 
admit an exterior complement that would solve the discrete counterpart of system (2.1a) and meet boundary 
condition (2.1b) in some approximate sense. 

We construct a second-order accurate discretization of system (2.1a) on the auxiliary Cartesian grid; 
a detailed description of the scheme can be found in Ref. [16]. The DPM will provide us with the com- 
plete boundary classification of all those and only those exterior grid vector-functions that solve the discrete 
counterpart of ( 2 . 1 a) outside the computational domain and meet boundary condition ( 2 . 1 b) (in the sense 
described below). The foregoing boundary classification will be obtained as an image of a special projection 
operator, which can be considered a discrete analogue to Calderon’s pseudodifferential boundary projec- 
tion [18, 19]. The projection operators act on the grid functions defined as boundary traces of the solution. 
In actual computations, the boundary conditions are set as follows. Every time we need to update the ghost 
cells we take an appropriate set of data from inside T (see below), project it onto the subspace in the entire 
space of boundary data that admits the correct exterior complement, and obtain the ghost cell values by 
calculating this complement on . 

The implementation of the DPM-based ABC’s starts with splitting the nodes of the auxiliary Cartesian 
grid into two distinct groups: those that are inside T and those that are outside T. Applying the stencil of 
the scheme for (2.1a) to each node of both groups, we consider the intersection of the grid sets swept by the 
stencil. This intersection is called the grid boundary 7; it is a multi-layered fringe of nodes of the auxiliary 
Cartesian grid located near and straddling the continuous boundary I\ 

For any function u on the Cartesian grid we define its trace Tr 7 u on 7 as merely a restriction. For 
any grid function u 7 specified on 7 we introduce the generalized potential Pw 7 with the density tz 7 ; the 
generalized potential is defined on the auxiliary Cartesian grid on 7 and outside it. The generalized potential 
is obtained as a solution of the special auxiliary problem (AP); solution of the AP replaces and extends the 
operation of convolution with the fundamental solution in classical potential theory. The AP is driven by 
the right-hand side that depends on n 7 , the formal construction of this right-hand side is the same in two- 
and three-dimensional cases, sec Refs. [7, 11, 12, 16] for details. The AP is formulated on a special finite 
auxiliary domain and the boundary conditions for the AP are chosen so that they approximate boundary 
condition (2.1b). The auxiliary domain is a Cartesian parallelepiped (i.e., aligned with the coordinate 
directions) that fully contains IV To make sure (2.1b) is taken into account properly, we specify periodicity 
boundary conditions in the cross- stream directions y and z. The periods are chosen sufficiently large to 
guarantee that the periodic solution considered on a finite fixed neighborhood of T and T 1 approximate well 
the theoretical non-periodic solution; the latter can be thought of as a limit when the periods in y and z 
approach infinity. The approximation of a non-periodic solution by the periodic one on a fixed subinterval 
as the period increases is discussed in our work [7, 11, 12]. 

Once the problem has been “periodized” in the cross-stream directions, we can separate the variables and 
then use a semi-analytic approach for the streamwise direction x. To do that, we apply the discrete Fourier 
transform in y and z to the finite- difference counterpart of (2.1a) and obtain a family of one- dimensional 
difference equations: 


5 



( 2 . 2 ) 


+ ^k^m—l,k ~ fm -1/2, ki 


Tfi 1, . , . , M, k — (fey, f’z)> 

ky 0, . . . , */y, fc 2 = 0, . . . , J z , 

where A k and B k are the 5 x 5 matrices and M + 1 , 2J y + 1 , and J z + 1 are the numbers of grid nodes in 
the x, y, and 2 directions, respectively (symmetry with respect to 2 = 0 is taken into account, as well as the 
fact that u and / are real- valued). Boundary conditions in the direction x are specified separately for each 
pair of wavenumbers k : 


(2.3a) 


S (fc) UqJc — 0, 


(2.3b) S+(k)ii M , k = 0, 

where 

(2- 4a ) S~(k) = n 


( 2 - 4b ) s+(k)= n (Qfc-^(fc)/), 

Qk = A k l B k , and Ms(fc) are the eigenvalues of Q k . The semi-analytic boundary conditions (2.3a) and 
(2.3b) (the eigenvalues for (2.4) are calculated numerically) explicitly prohibit growing modes of the solution 
in the left and right directions, respectively. 

In our work [12, 16], we have also discussed the possibility of replacing the Fourier transforms by non- 
unitary transforms. The latter may be needed when the grid in y and/or z is stretched (which provides for 
a drastic cost reduction) and consequently, the corresponding eigenfunctions form a skew basis. 

The foregoing AP allows us to calculate the generalized difference potential P u 7 for any grid density 
u y specified on 7. The composition of the operators Tr., and P, P 7 = Tr 7 P, is a projection, P/ = P 7 , 
and it is a discrete counterpart of Calderon’s boundary projection [18, 19] for system (2.1a). The image of 
this projection, ImP 7 , contains all those and only those u 7 ’s tin t are the traces of some exterior difference 
solution to (2.1a) that satisfies the boundary conditions of the AP — periodicity in y and 2 and boundary 
conditions (2.3) in x. The latter boundary conditions, in turn, approximate (2.1b). 

Having constructed the procedure for calculating the potent als and projections for the discrete version 
of (2.1a), we can now close the system inside the computationa domain, i.e., obtain the ABC’s. First, we 
take u and du/ ihi on 1', n is the normal, (these data arc availa )le from inside the computational domain) 
and, using interpolation liy along T and the first two terms of Taylor’s expansion (denoted 7t 7 ), obtain u 7 : 


(2.5) 


> \JL'y “ TC 



I 


Then, we need to calculate the potential P t> 7 for the density v 7 = P 7 u 7 and interpolate it to the nodes Fj : 
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( 2 . 6 ) 


u 


= Rr l Pv 1 = Firstly. 
Finally, the ABC’s are obtained in the operator form 


r 

where Tis composed of the operations of (2.5) and (2.6). Boundary condition (2.7) is applied every time we 
need to update the ghost cells when solving the interior problem (e.g., on every iteration). The implemen- 
tation of ABC’s (2.7) can either be direct or involve preliminary calculation of the matrix T. In the latter 
case, the runtime implementation of the ABC’s (2.7) is reduced to a matrix- vector multiplication. Numerical 
results for flows around the ONER A M6 wing obtained with the DPM-based boundary conditions (2.7) are 
summarized in work by Tsynkov and Vatsa [15] and Tsynkov [16, 17]. 

3. Application to Jet Flows. The major difference between the formulation of the previous section 
and the flow with jet exhaust is that in the vicinity of the jet we can no longer claim that flow perturbations 
against the free-stream background are small. Indeed, inside the propulsive jet the speed of the flow is 
typically much higher than the one in the surrounding area, moreover, other parameters, e.g., temperature, 
may also differ substantially. Therefore, the linearization of the flow against constant free-stream background 
that yields equations (2.1a), (2.1b) is, generally speaking, not valid in this case. 

However, in many applications (in particular, aerospace) one can clearly distinguish between those parts 
of the overall flow that contain jet(s) and the remaining areas. Therefore, the most comprehensive way 
to develop the far-field linearization in this situation will apparently be to use the appropriate asymptotic 
solutions for jets (see, e.g., Abramovich [26]) in the corresponding regions as a background. For flow regions 
outside the jet, it is always reasonable to assume that the foregoing linearization (2.1) will still be valid there. 

The particular setting that we will be studying hereafter is schematically shown in Figure 3.1. (The 
meaning of the two external grid surfaces is the same as T and T i in Figure 2.1.) It includes a three- 
dimensional slender body (symmetric with respect to the z = 0 plane but not axially symmetric, i.e., not a 
body of revolution) with sharp nose and boat-tail aft configuration; the rearmost plane surface of the body 
(not shown explicitly in Figure 3.1) is actually a location of the nozzle outlet; the outlet is rectangular in 
cross section. The exterior flow is subsonic with the free-stream Mach number Mo = 0.6 and zero angle of 
attack, the jet that is discharged from the outlet is supersonic, M 3 = 1.6, and confluent with the exterior 
flow. 

The specific shape of the body (see Figure 3.1) as well as parameters of the flow have been previously 
proposed for numerical study and analyzed in work by Compton [27]. In this original work [27], Compton had 
calculated external flow with the propulsive jet and also considered the interior portion of the flow, namely 
the flow in the actual nozzle located inside the afterbody (this nozzle flow obviously produces the jet). For 
our study, we have generated new grids and also simplified the overall formulation by eliminating the nozzle 
and specifying instead the uniform supersonic flow conditions at the nozzle outlet i.e., at the place where 
the jet starts. Compton’s goal [27] was to assess the performance of different turbulent models including 
their prediction capabilities for the flow inside the nozzle; our goal is to assess the performance of different 
external boundary conditions for the flows with jet exhaust. We, therefore, think that the aforementioned 
simplification is justified. 
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Block 1 



Fig. 3.1. Schematic geometric setup for the slender body with jet exhaust. 

Our typical grid consists of two blocks: block 1 of C-O type is for the exterior flow and block 2 of H-O 
type is for the jet portion (see Figure 3.1). Of course, this subdivision can only have an approximate meaning 
because the jet will obviously tend to spread while propagating downstream; basically, it means that the 
shear layer between the jet and coflow is located in the vicinity of the block interface. On this interface, the 
two grid blocks are point-matched, which is a requirement for TLNS3D. 

As has been mentioned, the exterior flow is subsonic and the jet is supersonic (other parameters of the 
flow will be pointed out later). The standard boundary condiiions in TLNS3D for this two- block jet flow 
case include one-dimensional characteristics for external inflow (block 1, upstream portion of the boundary), 
specification of the free-stream pressure with extrapolation of all other quantities for external outflow (block 
1, downstream portion of the boundary), extrapolation of all quantities for the jet downstream boundary 
(block 2) and specification of all quantities for the jet inflow boundary (block 2); the boundary conditions 
on the solid surface of the body are standard no-slip conditions. Extrapolation of all flow quantities at 
the jet outflow boundary is justified because as shown by numei ous simulations the core of the jet remains 
supersonic even at large distances downstream of the body, at least as far as 40 50 nozzle calibers away. 

The primary goal of this paper is to develop an alternative do the foregoing local boundary conditions 
for the jet flow case global ABC’s similar to those described in the previous section, and compare the 
performance of the two techniques. A direct implementation of lie ABC’s (2.7) will, however, encounter a 
major obstacle in this case: as has been mentioned, we cannot linearize against the free-stream background 
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in the jet region and therefore, cannot directly implement the operator T of (2.7) over the entire external 
boundary as this operator is obtained on the basis of the linear system (2.1). Of course, if we linearized 
the flow against constant free-stream background outside the jet and against some approximate asymptotic 
solution in the jet region (see Ref. [26]) and then used the corresponding linear system (unlike (2.1a), it will 
have variable coefficients) to construct the operator analogous to T of (2.7), then we could have applied 
the boundary conditions (2.7) straightforwardly as done in the previous work [15, 16] for flows with no 
jets. Computation of the new operator T in this framework will, in turn, require a different construction 
of the AP, certainly more elaborate (because of variable coefficients) and possibly more expensive than the 
one described in the previous section (see formulae (2.3), (2.4), (2.5)). One way of largely eliminating the 
difficulties associated with variable coefficients is apparently to take advantage of the supersonic nature of the 
jet and use marching-type algorithms in a subdomain of the new AP domain. Although this may make the 
whole foregoing program feasible, we consider its implementation as future work. In this paper we present 
the algorithm based on boundary conditions (2.7) with minimal alterations. 

As the ABC’s (2.7) obviously cannot be applied in the jet area, i.e., on that portion of the artificial 
boundary where the jet exits the domain, we need another procedure. The most natural choice will be 
to extrapolate all flow quantities downstream at the outflow boundary because the core of the jet remains 
supersonic even at large distances away from the nozzle outlet. Of course, we cannot actually predict where 
on the downstream boundary the flow actually becomes subsonic, i.e., where the extrapolation ceases to be 
applicable. However, we have observed that for the particular case under study we can extrapolate at least 
on the entire downstream boundary of the second grid block (see Figure 3.1). Thus, extrapolation of all flow 
quantities will be used henceforth as downstream boundary conditions for block 2. 

In the standard procedure, the downstream boundary conditions for grid block 1, i.e., on the rest of 
the outflow boundary, are based on the specification of free-stream pressure and extrapolation of all other 
quantities. Basically, these boundary conditions are relevant for subsonic outflow. In practice, some portion 
of the downstream boundary of block 1 may also be supersonic; in this case, however, the implementation of 
these pressure boundary conditions does not lead to noticeable errors because the streamwisc variations of 
pressure away from the nozzle are small (the jet is close to design, it is slightly overexpanded, see below) and 
therefore, specification of the free-stream pressure and extrapolation from the interior are both procedures 
with acceptable accuracy. 

To replace local boundary conditions on the outer boundary of block 1 (the region outside the jet) by the 
more accurate global ABC’s, we use relation (2.7). However, in formula (2.7) both the input and output are 
global, i.e., not only the operator T provides the ghost cell data along the entire boundary but also requires 
the data along the entire (penultimate) boundary as driving terms. By using extrapolation downstream in 
the jet core instead of using (2.7), we make sure that the possibly erroneous data from the global procedure 
are not used on this part of the boundary. However, as the global operator T is constructed on the basis of 
linearization (2.1), which is not valid in the jet area, plugging the actual flow quantities (including the jet 
profile) into the right-hand side of (2.7) may potentially generate errors along the entire outer boundary. 

On the other hand, it has been verified for model examples [14] and also seen for more complex cases 
that typically, closely located boundary nodes influence one another much stronger than the remote ones. 
This is a reasonable behavior from the standpoint of physics; in the structure of operators T it is reflected so 
that although the matrix is dense (non-locality) its near-diagonal terms are much larger than the off-diagonal 
ones (for systems as opposed to scalar equations, it will be a similar block- wise structure). The specific rate 
of decay for the off-diagonal terms can probably be obtained only for the most elementary formulations (e.g., 
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the Laplace equation with periodic boundary conditions). Hovever, although we cannot obtain analytical 
estimates for the kernels involved in operators T of (2.7), we can still make use of the actual (block- wise) 
off-diagonal decay in the numerical experiments. In practical terms, this implies that although substituting 
the jet profile into the right-hand side of (2.7) violates the small perturbations assumption , the associated 
error on the left-hand side of (2.7) will mostly be concentrated again in the jet area, where the results are 
not used for boundary conditions anyway as they are overridden by extrapolation. 

Thus, the actual combined DPM-based ABC’s that we employ for computation of the foregoing jet flow 
case are the following. For most of the outer boundary (except the near-jet area) we use formula (2.7) while 
substituting the actual flow profile in its right-hand side. For the jet core (outflow boundary for grid block 
2) we extrapolate all flow quantities downstream. For the small intermediate portion of the downstream 
boundary (near the jet core) we extrapolate all quantities except pressure, the latter is prescribed from 
its ffee-stream value. In fact, we have observed that within a certain range (5 to 30 cells of the fine grid 
described in Numerical Results), the actual location of where to switch from the pressure boundary condition 
to formula (2.7) does not exert much influence on either the final accuracy or multigrid convergence rate. In 
the next subsection, we provide an additional justification for the use of this procedure, 

3.1. Jet Outflow Boundary Conditions. To describe and explain the specific boundary treatment 
in the vicinity of the jet exit through the boundary, we start by considering the model problem below, 
disregarding for a moment the connection to the global boundary procedure described above. 

An model problem describing the error due to inaccurate outflow pressure data for the steady Euler 
equations is, 


(3.1) 


Ae x -|- Be y -1- Ce z =0, j < 0, 

p = 9{y,z), z = o, 

— oc < y, z < oo, 


where e — (p, u, v, w,p) T denotes the error and A , D , and C are constant matrices. We assume that the 
boundary data has compact support outside a small portion of the boundary, i.e. 

( 3 - 2 ) 9(y,z) = o, \y,z\>d, 

We also assume that the base flow is subsonic and moves to the light. The problem (3.1), (3. 2) is a model for 
the error in an approximate solution with correct outflow boundary data given on |y,z| > 5 and erroneous 
one’s on |y, z\ < 6. 

The relation of the model problem (3.1), (3. 2) to the specific Dutflow problem in this paper can briefly be 
described as follows. The global boundary procedure far away from the jet and the extrapolation procedure, 
sec Refs. [28, 29], in the supersonic part of the jet lead to very s nail errors, i.e. \g\ % 0. In an intermediate 
domain between the supersonic part of the jet and the part whe'e the global boundary conditions are used, 
pressure with erroneous data is specified, i.e. \g\ « 0(1) in that part of the domain. 

Note that for problems with boundary conditions in the (or streamwise) direction it makes little 
difference if one consider the inviscid Euler equations instead of tl e viscous thin layer Navier-Stokes equations 
since the number and nature of the boundary conditions require< i in the x direction are the same for the two 
sets of equations. 

Let Q" = ([”(£ + !)»—£] x R n_1 ) where n is the number of spatial dimensions, sec Figure 3.2. The 
following theorem describes the error distribution in the halfspa ;e x < 0. 
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Theorem 1. The error e in (3.1), (3.2) satisfies: 


(3.3) 


IHImos) ^ c . 


1 -H a 


\ \ f ( l -<? 2 


|Ap„ 


' <5 ' 


n — 1 


where C is a constant and |Ap max | is the maximal magnitude of the pressure error in \y>z\ < 8, Q = 
u 2 + v 2 + w 2 for n — 3 and Q 2 = u 2 - f- F 2 for n = 2. 



Fig. 3.2. The domain where L is t/ie distance /ram 17^ to the boundary. L » 1. 


Theorem 1 means that by measuring the error in a local L^-norm on the fixed domain £1*1, an error decay 
can be observed. Note that if the error was measured by computing the L 2 norm in the whole computational 
domain , no error decay could be observed. The proof of theorem 1 involves a straightforward application of 
the theory of Ref. [30]. Numerical experiments that verify the decay rate (3.3) can be found in Ref. [31]. 

For our specific outflow problem with erroneous data given on the intermediate domain between the 
supersonic part of the jet and the part where global boundary conditions are used, theorem 1 means that the 
error decays with the rate 6 2 / L away from the outflow boundary. Furthermore, in our specific flow problem 
we have a slightly overexpanded jet which means that the maximum pressure error |Ap max | in (3.3) is rather 
small. 

3.2. Effective Reynolds Number. To calculate the operator T of (2.7), we are solving the AP for 
system (2.1). This system is obtained by linearization of the original thin-layer equations. However, as 
the actual flows that we are studying are turbulent, to integrate the thin-layer equations numerically one 
complements them with turbulence models inside the computational domain. These models may be complex 
and require solving some additional differential equations (see next section). 

For the simplified linearized far- field representation, we do not use these accurate and sophisticated 
turbulence models. However, we still need to account for the corresponding turbulent mixing and dissipation, 
at least in an approximate way. In the previous work [9] Tsynkov, et al. have used the concept of effective 
turbulent viscosity for the far field and calculated the effective turbulent Reynolds number using the fact 
that the laminar and turbulent plane wakes behind the body have the same asymptotic behavior. [32] 
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The asymptotic behavior of laminar and turbulent circular jets is also known to be the same [26, 32]. It 
involves a linear increase in width and a decrease in center-line velocity inverse proportional to the distance 
from the source. The virtual kinematic viscosity (incompressib e case) can be considered constant through 
the entire jet region. Although we do not use boundary conditions (2.7) in the core of the jet, the outer 
portions of the shear layer region are still covered by the global procedure, therefore we need to provide the 
effective value of 1/Re in equations (2.1). 

The jet that we are studying is rectangular in its initial cross section (see next section for particular 
details); however, its shape will approach circular further away of the outlet. Therefore, we will use the 
results obtained for circular jets to find an approximate value for the effective Reynolds number. First, we 
notice that the universal velocity profiles in a cross section of an incompressible submerged jet (i.e., the 
jet that propagates through a medium at rest) are the same as those obtained for the excess velocity of 
the jet propagating in a coflow. [26] Moreover, many experimental observations corroborate [26] that the 
same universal profiles remain valid for a compressible supersonic jet spreading through either a stationary 
or moving medium. Of course, while the profiles are universal, the actual spreading rate for the jet will 
differ for different cases. Second, for the particular case under study (the ratio of stagnation temperatures 
is T* /Tq = 0.936; the design pressure ratio is Pj/Po\ design = 4.25 at Mj = 1.6 whereas the actual pressure 
ratio is p*/p 0 = 4.00, the jet is slightly overexpanded), the initial value of the compressibility parameter [26] 
is p = pq/pj = 1.41 and the initial velocity ratio is m — u Q /u 7 = ^Tq/TjMo/Mj = 0.459. According to 
Ref. [26], these values are within the range (0 < m < 0.6, 0.3 < p < 1.43), for which the correction due to 
compressibility for the spreading rate b of the jet can be taken into account by calculating it as 

(3.4a) 

instead of the old expression 


J comp 


lfp 1 - m 
2 1 + pm 


( 34b ) bine = Cl 1 --” 1 , 

1 + m 

which is relevant for the incompressible flow; c in formulae (3.4) is a constant and x is the distance from the 
source. 

According to the measurements referenced by Schlichting [32], for a submerged incompressible jet 6,/ 2 = 
0.0848x, where 6 j/ 2 is half width of the jet at half depth. Substituting this into the solution for laminar 
jet [26, 32]: 

ftl/2 

X 

one obtains the virtual kinematic viscosity [32] : 


5.27 


yfK' 


(3.5) ur = 0.0161 y/E, 

here K is the total kinematic momentum flux. Since the velo ;ity profiles are universal, for the jet with 
coflow we only need to multiply the spreading rate by (1 - m)/( 1 -f m) according to formula (3.4b) and for 
the compressibility correction we use (3.4a), which altogether yields: 
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(3.6) 


1ST — 0.00636 y/K. 


As has been mentioned, the boundary condition that we specify for the jet inflow is a uniform supersonic 
profile across the entire nozzle outlet. Therefore, the quantity K can be obtained by multiplying the square 
of the excess velocity (relative velocity of the jet with respect to the velocity of coflow) by the area of the 
outlet a, K = ( Uj — uo) 2 &- Then, the effective turbulent Reynolds number is calculated as Rer = UL/is t , 
where U is the characteristic speed and L is the characteristic length. For the particular setting under study, 
it is reasonable to assume that U = \ Uj — uq| and L = \fo. Consequently, from (3.6) we conclude that 


(3.7) Re r = 0.00636 -1 » 157. 

In our computations, the actual value of Re for system (2.1) was taken from (3.7). 

4. Numerical Results. The particular geometry of the body shown in Figure 3.1 is the following: 
rectangular cross section y x z = 6.2 x 6.8 with rounded edges; sharp nose and boat-tail afterbody; total 
length in the x direction is 63; rectangular nozzle outlet y x z — 2.62 x 5.12; full description can be found in 
the work by Compton [27]. 

The geometry and the flow are symmetric with respect to the plane z = 0 (zero angle of attack). For 
our computations we have used three different domains with successively reduced dimensions, see Figure 4.1; 
domain I (or large domain) with the diameter of about 30 calibers of the body was used for calculating the 
reference solutions, domain II is 0.36 or about 1/3 of the size of domain I in each direction and domain III 
is 0.22 or about 1/5 of the size of domain I in each direction. 
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As has been mentioned, to integrate the thin- layer equation: on the curvilinear grid shown on Figure 3.1 
we use the code TLNS3D by Vatsa, et al. [23, 24] This is a central-difference code with five stage explicit pseudo- 
time Rungc-Kutta relaxation used for obtaining steady-state solutions. The code employs local Courant step, 
semi-implicit residual smoothing, and multigrid for accelerating the convergence. In our computations, we 
used either three or two nested grid levels with V cycles (depending on the grid dimension); this multi-level 
V-cycle algorithm is, in fact, a final stage of the full multigrid (FMG) procedure. In addition, to improve 
the convergence to steady state, the solver is preconditioned according to the methodology of Ref. [33], 

The DPM-based ABC’s arc implemented only on the finest grid for the V-cycle in the final FMG stage; 
the boundary data for coarser levels are provided by the coarsening procedure. Moreover, even on this finest 
grid we implement the DPM-based ABC’s only on the first and the last Runge-Kutta stages, which has been 
found (15, 16] to make very little difference compared to the implementation on all five stages; the boundary 
data for the three intermediate stages are provided from the DPM-based ABC’s on the first stage. 

To account for the turbulent phenomena, the solver is also supplemented with Menter’s two-equation 
turbulence model [34], The actual molecular Reynolds number based on unit length is Re = 321000, Prandtl 
number is Pr = 0.72, specific ratio is « = 1.4. 

We have used several different grids to calculate the jet flew; in all cases we kept the normal spacing 
near the solid surface the same: « 3 • 10 -4 . All grids are stretched, the cell size increases away of the body in 
geometric progression. The dimension of the C-0 grid block 1 for domain I was ixjxk= 385 x 77 x 33 (i is 
the streamwise C-type coordinate, j is the radial coordinate, and k is the circumferential cross-stream O-typc 
coordinate, quarter circle). The dimension of the H-0 grid block 2 for domain I was i x j x k = 81 x 77 x 65 
(* is streamwise, j is radial, and k covers half circle). We will further refer to this grid as fine. On the 
fine grid, we have calculated two reference solutions, one with standard ABC’s and another with global 
ABC’s. As the artificial boundary for domain I is located sufficiently far away of the body, the diff erence 
between the corresponding results is negligible. On Figures 4.2, we show convergence histories for this case: 
residual of the continuity equation is plotted vs. work units on Figure 4.2a and drag coefficient is plotted vs. 
work units on Figure 4.2b. (One work unit is the cost of advancing one time step on the finest grid.) 


Convergence histories for the jet flow 
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Fig. 4.2a. Convergence histones for the residual of the continuity equation, domain I, fine gnd. 


14 



Drag convergence for the jet flow 
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Fig. 4.2b. Drag convergence, domain I, fine grid. 


^From Figure 4.2 we conclude that multigrid convergence rates arc the same for local and global ABC’s 
on domain I. Moreover, the values of total drag coefficient per unit area Co are summarized for this case in 
the right column of Table 4.1. They differ by about one third of one per cent, which corroborates that the 
type of external boundary conditions has little effect on the solution itself, as well as multigrid convergence 
history, for large computational domains. 


Table 4.1 

Total drag coefficient per unit area Co- 


Domain 

hi 

II 

I 

Grid 

fine 

coarse 

fine 

fine 

Local 

— 

2.77 ±.03 

2.74 ± .04 

2.506 

Global 

2.365 

2.495 

2.484 

2.497 


For domain II, we have computed the flow on two grids with different dimensions. The first grid has the 
same number of nodes as the one used in domain I; it was, in fact, constructed by scaling down the original 
large grid by a factor of 0.36 in each direction. We will also refer to it as fine grid. As shown in Table 4.1, 
the coefficient Co obtained on this grid with global ABC’s differs by less than one per cent from its reference 
value, whereas the accuracy provided by local ABC’s is not nearly as good, about 9% discrepancy; moreover, 
because of the poor convergence (see Figures 4.3) the value of Co for local ABC’s is given with the error 
bands indicated. 

The much smaller size of domain II compared to domain I actually suggests that on domain II one can 
successfully compute the solution on a grid with fewer nodes. Therefore, the second grid that we have used for 
domain II had one half of the original dimension in two out of three directions, block 1 ixjxk = 193 x 39 x 33 
and block 2ixjx/c = 41x39x 65, this grid will be referred to as coarse . Again, as follows from Table 4.1, 
global ABC’s provide for an accurate solution whereas the accuracy of local ABC’s is not sufficient and the 
convergence is slow (or even non-existent). Convergence histories for domain II are presented on Figures 4.3. 

Since the node count for the coarse grid is only 1/4 of the node count for the fine grid, the convergence 
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vs. work for the coarse grid is about four times faster (see Figures 4.3), although convergence rates measured 
vs. number of multigrid cycles will be approximately the same for both grids. Note that because of the 
particular grid dimensions (the issue of divisibility by 2) we have used three nested multigrid levels on the 
fine grid and two levels on the coarse grid. One can clearly see from Figures 4.3 that the DPM-based ABC’s 
provide for a noticeably higher multigrid convergence rate than the standard local ABC’s do. Moreover, it 
is, in fact, hard to conclude from Figures 4.3 whether or not the algorithm with local ABC’s converges. If it 
does, the resulting C D will be about 10% off its reference value. 


Convergence histories for the jet flow 
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Fig. 4.3a. Convergence histones for the residual of the continuity equation, domain II, fine and coarse grids. 


Drag convergence for the jet flow 



Fig. 4.3b. Drag convergence, domain II, fin ? and coarse grids. 

On domain III, the computations were performed on the fin< grid, which again was obtained by scaling 
down the grid from domain I (a factor of 0.22 in each direction}. The algorithm with local ABC’s for this 
domain/grid failed to converge, whereas the algorithm with global ABC’s converged with the same rate 
as before. However, the actual computed Co is about 5% off its reference value (see Table 4.1). T his 
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can apparently be attributed to the fact that the assumption of linearity (small perturbations) outside the 
computational domain is violated for such a small domain size. Convergence histories for domain III arc 
presented on Figures 4.4. 


Convergence histories for the jet flow 
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Fig. 4.4a. Convergence histories for the residual of the continuity equation , domain III, fine grid. 


Drag convergence for the jet flow 
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Fig. 4.4b. Drag convergence, domain III, fine grid. 

Computations on a coarse grid for domain III were not performed because we did not expect to recover 
the accurate value of Cp- However, the fact that the algorithm with global ABC’s converges on domain III 
corroborates the high robustness of this procedure. 

All computations described in this section were conducted on Cray Research computers, J90 and C90 
series. Computational overhead due to the use of global ABC’s is about 15% for the particular fine grid 
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referenced before. This overhead is determined mostly by dox lain geometry and typically does not scale 
linearly with the dimension of interior grid. For the aforementi< 'ned coarse grid the overhead reaches 30%. 

5. Conclusions. We have constructed and implemented global ABC’s for calculating external flows 
with jet exhaust. The ABC’s combine extrapolation of all flow quantities downstream in the supersonic 
core of the jet and nonlocal DPM- based treatment for the remaining portion of outer boundary. The 
overhead associated with implementation of the new technique is is compensated for by the reduced grid 
dimension on small domains and higher convergence rate. In the series of computations performed, the 
DPM-based algorithm have consistently demonstrated better accuracy, faster multigrid convergence, and 
higher robustness compared to the standard local methodology. 

6. Acknowledgment. We arc most grateful to Dr. E. B. Parlette of Vigyan, Inc., for his valuable help 
in grid generation. 
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